      subroutine current
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M, ONLY:
      use cophys_com_M
      use ctemp_com_M, ONLY: fmf, fef, fvf, fbxf, fbyf, fbzf, twx
      use blcom_com_M
      use objects_com_M, ONLY: chi
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
!
!     A subroutine that uses curl(B)=4pi/c J
!     to get J from B
!
!     called by vinit.f and by celest3d
!

      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer , dimension(8) :: iv
      integer, dimension(4) :: lioinput, lioprint
      integer :: irg
      real(double), dimension(1) :: ixi1, ixi2, ieta1, ieta2
      real(double) :: ixi4, ieta4, ixi5, ieta5
      real(double), dimension(8) :: wght
      real(double), dimension(100) :: fpxf, fpyf, fpzf, fpxft, fpyft, fpzft, &
         ucm, vcm, wcm, cm
      real(double), dimension(3,20) :: wsin, wcos
      real(double), dimension(3,12,20) :: tsi, rot
      real(double) :: dummy, factor, wsp
      logical :: expnd, nomore
      character :: name*80
!-----------------------------------------------
!
!	Compute J0 =curl(B) c/4pi
!

      call curlv (nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, c6x&
         , c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z&
         , c5z, c6z, c7z, c8z, bxn, byn, bzn, jx0, jy0, jz0)
!
!     impose toroidal boundary conditions on j0
!

      dummy = 0.0D0

      call torusbc (ibp1 + 1, jbp1 + 1, kbp1 + 1, cdlt, sdlt, dummy, dz, jx0, &
         jy0, jz0, periodicx)
!
      call axisgrad (ibp1, jbp1, iwid, jwid, kwid, jx0, jy0, jz0)
!
      factor = fourpi/clite
!
      wsp = 1.
      call torusbc_scal (2, ibp2, 2, jbp2, 2, kbp2, iwid, jwid, kwid, wsp, chi&
         , periodicx)
!
      do k = 1, kbp2
         do j = 1, jbp2
            i = 1
            if (ibp2 > 0) then
!
!
!
               jx0((k-1)*kwid+(j-1)*jwid+1:(ibp2-1)*iwid+(k-1)*kwid+(j-1)*jwid+&
                  1:iwid) = jx0((k-1)*kwid+(j-1)*jwid+1:(ibp2-1)*iwid+(k-1)*&
                  kwid+(j-1)*jwid+1:iwid)*factor
               jy0((k-1)*kwid+(j-1)*jwid+1:(ibp2-1)*iwid+(k-1)*kwid+(j-1)*jwid+&
                  1:iwid) = jy0((k-1)*kwid+(j-1)*jwid+1:(ibp2-1)*iwid+(k-1)*&
                  kwid+(j-1)*jwid+1:iwid)*factor
               jz0((k-1)*kwid+(j-1)*jwid+1:(ibp2-1)*iwid+(k-1)*kwid+(j-1)*jwid+&
                  1:iwid) = jz0((k-1)*kwid+(j-1)*jwid+1:(ibp2-1)*iwid+(k-1)*&
                  kwid+(j-1)*jwid+1:iwid)*factor
!
               i = ibp2 + 1
               ijk = (ibp2 - 1)*iwid + (k - 1)*kwid + (j - 1)*jwid + 1
            endif
         end do
      end do
!
!	Add an extgernal current chosen to zero out the initial current of the plasma species
!	Needed to run eletrostatic-like cases such as the two stream instability
!
	if(external_current) then

	   do irg=1,nrg
	      jx0=jx0+rhr(irg)*uvi(irg)*qom(irg)/abs(qom(irg))*vvol
	      jy0=jy0+rhr(irg)*vvi(irg)*qom(irg)/abs(qom(irg))*vvol
	      jz0=jz0+rhr(irg)*wvi(irg)*qom(irg)/abs(qom(irg))*vvol
	   enddo

	end if
      return
      end subroutine current
